#code to replicate analyses in Hainmueller (2012) adapted from replication materials 
#for https://doi.org/10.1093/pan/mpr025 
rm(list=ls())
library(Matching)
library(survey)
library(lattice)
library(ebal)
library(rgenoud)
library(haven)
library(tidyverse)
library(grid)
library(gridExtra)
library(cowplot)

dat <- read_dta("data/output/iraqmerged.dta")

#remove Kurdish respondents and set covars
dat <- dat %>%
  filter(kurd==0) %>%
  dplyr::select(treat,natident,age, gender, coledu, 
                pctshia, sqsect, clanperthou)

dat <- dat[order(dat$treat),]

colnames(dat) <-gsub("age","Age",colnames(dat))
colnames(dat) <-gsub("gender","Gender (1=female)",colnames(dat))
colnames(dat) <-gsub("coledu","College educ.",colnames(dat))
colnames(dat) <-gsub("pctshia","Shi'i %",colnames(dat))
colnames(dat) <-gsub("sqsect","Prehist. vio.",colnames(dat))
colnames(dat) <-gsub("clanperthou","Tribal ties",colnames(dat))

dat <- dat%>%
  na.omit()

covarsname <- colnames(dat)    
X <- dat[,covarsname]
W <- as.integer(X$treat)
Y <- as.integer(X$natident)
X <- X[,-c(1:2)]
cvars <- colnames(X)
X <- as.matrix(X)

# balance before matching
bout.nm  <- MatchBalance(W~X,match.out = NULL,ks=FALSE)
bal.nm   <- baltest.collect(matchbal.out=bout.nm ,var.names=colnames(X),after=FALSE)
round(bal.nm,2)

# Maha Dist Matching
mout.maha  <-  Match(Y,W,X,BiasAdjust=F,estimand="ATT",M=1)
summary(mout.maha)
bout.maha  <- MatchBalance(W~X,match.out = mout.maha,ks=FALSE)
bal.maha   <- baltest.collect(matchbal.out=bout.maha ,var.names=colnames(X),after=TRUE)
round(bal.maha,2)

# Genetic Matching
g.weights <- GenMatch(Tr=W, X=X, BalanceMatrix=X, estimand="ATT", M=1,print.level=0)
mout.gm   <-  Match(Y,W,X,BiasAdjust=F,Weight.matrix=g.weights,estimand="ATT",M=1)
summary(mout.gm)
bout.gm <- MatchBalance(W~X,match.out = mout.gm,print.level=0,ks=FALSE)
bal.gm <- baltest.collect(matchbal.out=bout.gm,var.names=colnames(X),after=TRUE)
round(bal.gm,2)

## Logistic PS weighting + matching
PS  <- glm(W~X,family=binomial(link=logit))$fitted
PSM <- PS
PS  <- PS[W==0]
PS  <- PS/(1-PS)

# PS Matching 
mout.psm <-  Match(Y,W,X=PSM,BiasAdjust=F,estimand="ATT",M=1)
summary(mout.psm)
bout.psm <- MatchBalance(W~X,match.out = mout.psm,print.level=0,ks=FALSE)
bal.psm <- baltest.collect(matchbal.out=bout.psm,var.names=colnames(X),after=TRUE)
round(bal.psm,2)

# PS Weighting
bout.psw <- MatchBalance(W~X,weights=c(PS,W[W==1]),ks=FALSE)
bal.psw <- baltest.collect(matchbal.out=bout.psw,var.names=colnames(X),after=FALSE)
round(bal.psw,2)

# Entropy Balancing
out.eb <- ebalance(
  Treatment=W,
  X=X
)

bout.eb <- MatchBalance(W~X,weights=c(out.eb$w,W[W==1]),ks=FALSE)
bal.eb  <- baltest.collect(matchbal.out=bout.eb,var.names=colnames(X),after=F)
round(bal.eb,2)

# # Entropy Balancing (with trimmed weights)
# out.ebtr <- ebalance.trim(          
#   ebalanceobj=out.eb
# )
# 
# bout.ebtr <- MatchBalance(W~X,weights=c(out.ebtr$w,W[W==1]),ks=FALSE)
# bal.ebtr  <- baltest.collect(matchbal.out=bout.eb,var.names=colnames(X),after=FALSE)
# round(bal.ebtr,2)
# 
# dat$w   <-  c(out.eb$w,W[W==1])
# dat$wtr <-  c(out.ebtr$w,W[W==1])

# Effect Estimates
# des <- svydesign(id=~1,weights=~w, data= dat)
# mod1 <- svyglm(vote_l_97 ~ tolabor, design = des)
# summary(mod1)
# 
# des <- svydesign(id=~1,weights=~wtr, data= dat)
# mod2 <- svyglm(vote_l_97 ~ tolabor, design = des)
# summary(mod2)

# Balance Plots
d  <- rbind(bal.nm,bal.maha,bal.gm,bal.psm,bal.psw,bal.eb)
dn <- rownames(d)
d  <- data.frame(d)
d$vname <- factor(dn,levels = unique(dn)[length(unique(dn)):1], labels = unique(dn)[length(unique(dn)):1])
d$gr <- rep(c("Unadjusted","MahaDist Matching","Genetic Matching","PS Matching","PS Weighting","Entropy Balancing"),each=length(unique(dn)))
d$gr <- factor(d$gr,levels=unique(d$gr)[1:length(unique(d$gr))],labels=unique(d$gr)[1:length(unique(d$gr))])

mypal<-c("black",rep("darkgrey",4),"black")[6:1]
Cex <- Cex2 <- 1

# plot with SDs
bplot <- function(x,y,...)
{
  panel.abline(v=0, lwd = 1 , lty="solid")
  panel.abline(v=c(-.1,.1), lwd = 2 , lty="dotted")
  panel.abline(h=c(1:nrow(d)), lwd = 1 , lty="dashed", col="gray95")
  panel.xyplot(x,y,...)
}

png("plots/compmatch1.png", 
    width = 200, height = 150, units='mm', res = 300)
print(
  
  xyplot(vname~round(sdiff.pooled/100,2),data=d,groups=gr,xlim=c(-.6,.6),
         panel = bplot
         ,par.settings = list(superpose.symbol = list(pch = c(15,19,17,18,15,1)[6:1],col=mypal,cex=1))
         ,xlab=list("standardized difference in means",cex=Cex),ylab="",auto.key = list(columns=3),scales=list(y=list(cex=Cex),x=list(cex=Cex2,at=c(-.5,-.1,0,.1,.5),labels=c("-.5","-.1","0",".1",".5")))
  )
)
dev.off()

png("plots/compmatch2.png", 
    width = 200, height = 150, units='mm', res = 300)
# plot with p-values
print(
  xyplot(vname~round(T.pval,2),data=d,groups=gr,xlim=c(-.05,1.05),
         panel = bplot
         ,par.settings = list(superpose.symbol = list(pch = c(15,19,17,18,15,1)[6:1],col=mypal,cex=1))
         ,xlab=list("p-value: difference of means test",cex=Cex),ylab="",auto.key = list(columns=3),scales=list(y=list(cex=Cex),x=list(cex=Cex2,at=c(0,.1,.2,.5,1),labels=c("0",".1",".2",".5","1")))
  )
)
dev.off()

g <- as_grob(print(
  xyplot(vname~round(sdiff.pooled/100,2),data=d,groups=gr,xlim=c(-.6,.6),
         panel = bplot
         ,par.settings = list(superpose.symbol = list(pch = c(15,19,17,18,15,1)[6:1],col=mypal,cex=1))
         ,main = "standardized difference in means",xlab=list("",cex=Cex),ylab="",auto.key = F,
         scales=list(y=list(cex=Cex),x=list(cex=Cex2,at=c(-.5,-.1,0,.1,.5),labels=c("-.5","-.1","0",".1",".5")))
  )
))

g1 <- as_grob(print(
  xyplot(vname~round(T.pval,2),data=d,groups=gr,xlim=c(-.05,1.05),
         panel = bplot
         ,par.settings = list(superpose.symbol = list(pch = c(15,19,17,18,15,1)[6:1],col=mypal,cex=1))
         ,main = "p-value: difference of means test", xlab=list("",cex=Cex),ylab="",auto.key = list(space="bottom",columns=3, border=F),
         scales=list(y=list(cex=Cex),x=list(cex=Cex2,at=c(0,.1,.2,.5,1),labels=c("0",".1",".2",".5","1")))
  )
))

plot_grid(g1, g, labels = "AUTO", rel_heights = c(1,1.1), nrow=2)

###Figure B.3###
png("plots/compmatchb.png", 
    width = 240, height = 160, units='mm', res = 300)
plot_grid(g, g1, labels = "AUTO", rel_heights = c(1,1.1), nrow=2)
dev.off()
